#library(ggplot2)
#theme_set(theme_classic() +
#		theme(text=element_text('times')) +
#		theme(strip.background=element_blank()) +
#		theme(strip.text.x=element_text(face='bold',size=12)))
#library(grid)  # unit

# factor(round(runif(10,0,3)), levels=0:3, labels=c('T','A1','A2','D'))

coutinho <- function(tmax, L=350, showimage=FALSE, benchmark=FALSE)
{
	Phiv <- 0.05
	tau <- 4
	Pinf <- 4e-5 #1e-5
	R <- 4
	Prepl <- 0.99
	
	rule1 <- function(nA1, nA2, t)  # T
		ifelse(nA1 >= 1 | nA2 >= R, 1, 0)
	rule2 <- function(nA1, nA2, t)  # A1
		ifelse(t >= tau, 2, 1)
	rule3 <- function(nA1, nA2, t)  # A2
		3
	rule4 <- function(nA1, nA2, t)  # D
		ifelse(rbinom(1,1,Prepl), ifelse(rbinom(1,1,Pinf), 1, 0), 2)

	graphic <- function(ca) {
		print(t)
		noaxis <- theme(axis.line=element_blank(),
			axis.text.x=element_blank(),
			axis.text.y=element_blank(),
			axis.ticks=element_blank(),
			axis.title.x=element_blank(),
			axis.title.y=element_blank(),
			legend.position="none",
			panel.background=element_blank(),
			panel.border=element_blank(),
			panel.grid.major=element_blank(),
			panel.grid.minor=element_blank(),
			plot.background=element_blank(),
			plot.margin=unit(c(-0.04,-0.05,-0.08,-0.06),'npc'))  #(top,right,bottom,left)
		l <- 20  # because of memory constrains, lets cut size
		df <- expand.grid(x=1:l, y=1:l)
		df$z <- factor(apply(df, 1, function(d) ca[d[1],d[2]]), levels=0:3)
		p <- ggplot(df, aes(x,y,fill=z)) +
			geom_tile() +
			scale_fill_manual(values=c('green','yellow','red','gray')) +
			noaxis + labs(x=NULL, y=NULL)
		#ggsave(paste0('tick',t,'.png'), p, png, width=l*2, height=l*2, limitsize=FALSE)
		print(p)
		p
	}
	
	par(mar=c(0, 0, 0, 0), xaxs='i', yaxs='i')

	graphic2 <- function(ca) {
		# ggplot2 too slow -- using base
		image(ca, axes=FALSE, main="", xlab="", ylab="", zlim=c(0,3), col=c('green','yellow','red','gray'), useRaster=TRUE)
	}

	# prebuild data frame results
	# http://stackoverflow.com/questions/19697700/how-to-speed-up-rbind
	ca <- matrix(0, L,L)
	ca[sample(1:(L*L), L*L*Phiv)] <- 1
	nA1 <- matrix(0, L,L)
	nA2 <- matrix(0, L,L)
	time <- matrix(0, L,L)
	fn <- list(rule1,rule2,rule3,rule4)

	temp <- rep(0,tmax+1)
	X <- data.frame(t=temp, T=temp, A=temp, D=temp)
	t <- 0
	if(benchmark)
		tic <- proc.time()
	while(t <= tmax) {
		#if(any(t == c(0:10,18:20,25,35,45,50,80,150,200)))
		if(showimage && t==52)
			graphic2(ca)
		X[t+1,] <- c(t, sum(ca==0), sum(ca==1)+sum(ca==2), sum(ca==3))

		# two phases so that the update of the headCA happens synchronously

		for(x in 0:(L-1))
			for(y in 0:(L-1))
				if(ca[x+1,y+1] == 0) {
					get <- function(dx, dy)
						ca[(x+dx)%%L+1,(y+dy)%%L+1]
					# neighborhood configuration from the paper
					N <- c(get(-1,0),get(-2,0),get(+1,0),get(+2,0),
						get(0,-1),get(0,-2),get(0,+1),get(0,+2))  # Cross r=2
					#N <- c(get(-1,-1),get(-1,0),get(-1,+1),get(0,-1),
					#	get(0,+1),get(+1,-1),get(+1,0),get(+1,+1))  # Moore r=1
					nA1[x+1,y+1] <- sum(N==1)
					nA2[x+1,y+1] <- sum(N==2)
				}

		for(x in 0:(L-1))
			for(y in 0:(L-1)) {
				n <- ca[x+1,y+1]
				m <- fn[[n+1]](nA1[x+1,y+1], nA2[x+1,y+1], time[x+1,y+1])
				if(m != n) {
					ca[x+1,y+1] <- m
					time[x+1,y+1] <- 0
				}
				else
					time[x+1,y+1] <- time[x+1,y+1]+1
			}
		#if(t > 14)
		#	readline(paste0("t=",t,". Press <return> to continue"))
		t <- t+1
	}
	if(benchmark)
		return(sum((proc.time()-tic)[c(1,4)]))
	X
}


draw <- function(X, filename)
{
	library(ggplot2)
	library(reshape2)
	library(plyr)
	theme_set(theme_classic() +
			theme(text=element_text('serif')) +
			theme(strip.background=element_blank()) +
			theme(strip.text.x=element_text(face='bold',size=12)))

	X <- subset(X, t<=350)
	Y <- melt(X, id.vars='t')
	Y$value <- Y$value/max(Y$value)

	Y[Y$t>12,'t'] <- Y[Y$t>12,'t'] + 120-12
	Y[Y$t<=12,'t'] <- Y[Y$t<=12,'t']*10

	p <- ggplot(Y, aes(x=t, y=value, color=variable, group=variable)) +
		geom_line(subset=.(t<=120)) +
		geom_line(subset=.(t>120), size=1) +
		geom_point(subset=.(t<=120), size=2) +
		geom_vline(xintercept=12*10, linetype=2, size=0.4) +
		scale_color_manual("Cell State", values=c('green','orange','gray'),labels=c('T','A1+A2','D')) +
		xlab('iterations') + ylab('Cell Density') +
		scale_x_continuous(breaks=c(seq(0,12,4)*10, seq(50,350,50)+120-12), labels=c(seq(0,12,4),seq(50,350,50))) +
		#scale_x_continuous(breaks=c(c(0,seq(3,scalet,3)*52, seq(52+52^2,300,52))), labels=c('0',paste(seq(3,scalet,3),'weeks'), '1 year', paste(2:3,'years'))) +
		theme(legend.position='none')
	ggsave(paste0(filename,'.pdf'), p, width=6, height=3)
}
